---
title: 比较 stats::procmp() 和 vegan::rda() 两种方法
author: Gao
date: '2018-11-02'
slug: comparison-of-two-pca-methods
categories:
  - R
tags:
  - PCA
  - statistics
---

<script src="/rmarkdown-libs/header-attrs/header-attrs.js"></script>


<p>上次提及 PCA 分析的方法有很多种。那不同方法之间的得到的结果会有差异吗？</p>
<p>最近采用 PCA 分析 RNA-seq 样本之间的差异，得到了下面的结果。</p>
<div id="生成示例数据" class="section level1">
<h1>生成示例数据</h1>
<p>生成一个含有 1000 个基因, 27 个样品的数据集.
这 27 个样品来自于 3 个基因型(WT, Mutant1, Mutant2),
3 种处理(CK, Trt1, Trt2), 共分为 <span class="math inline">\(3 * 3 = 9\)</span> 组,
每组 3 个重复, 合计 27 个样品.</p>
<pre class="r"><code>library(DESeq2)
library(dplyr)
library(ggplot2)
theme_set(theme_bw())

dds &lt;- makeExampleDESeqDataSet(n = 1000, m=27)
dds$condition &lt;- factor(rep(rep(c(&quot;Ctrl&quot;,&quot;Trt1&quot;,&quot;Trt2&quot;),3), 3))
dds$genotype &lt;- factor(rep(rep(c(&quot;WT&quot;,&quot;MuA&quot;,&quot;MuB&quot;), each=3),3))

# sample table
sample_table &lt;- colData(dds) %&gt;% as.data.frame() %&gt;%
  tibble::rownames_to_column(var=&quot;sample_id&quot;)

# 做 log2 变换
rld &lt;- rlog(dds, blind = F)</code></pre>
</div>
<div id="使用-statsprcomp-进行主成分分析" class="section level1">
<h1>使用 <code>stats::prcomp()</code> 进行主成分分析</h1>
<pre class="r"><code># 运行 PCA
pca &lt;- stats::prcomp(t(assay(rld)))

# 计算解释度
percent_var &lt;- pca$sdev^2/sum(pca$sdev^2)

# 绘图
df &lt;- data.frame(PC1=pca$x[,1], 
PC2=pca$x[,2], sample_id=colnames(rld)) %&gt;% 
  left_join(sample_table)

mapping &lt;- aes(PC1, PC2, color=condition, shape=genotype,label=sample_id)

p1 &lt;- ggplot(df,mapping) +
    geom_point(size=3) +
    ggrepel::geom_text_repel(show.legend = F) +
    xlab(paste0(&quot;PC1: &quot;, round(percent_var[1] * 100), &quot;% variance&quot;)) + 
    ylab(paste0(&quot;PC2: &quot;, round(percent_var[2] * 100), &quot;% variance&quot;)) +
    ggtitle(&quot;method: stats::prcomp()&quot;)
    
p1</code></pre>
<p><img src="/post/2018-11-02-comparison-of-two-pca-methods_files/figure-html/pca1-1.png" width="768" /></p>
</div>
<div id="使用-veganrda-做主成分分析" class="section level1">
<h1>使用 <code>vegan::rda()</code> 做主成分分析</h1>
<pre class="r"><code># rda() 分析
pca &lt;- vegan::rda(t(assay(rld))) 

# 计算解释度
percent_var &lt;- pca$CA$eig/pca$tot.chi  # rda() 的结果中信息比较完整
df &lt;- vegan::scores(pca)$sites %&gt;% 
  as.data.frame() %&gt;%
  tibble::rownames_to_column(var=&quot;sample_id&quot;) %&gt;%
  left_join(sample_table)

p2 &lt;- ggplot(df,mapping) +
    geom_point(size=3) +
    ggrepel::geom_text_repel(show.legend = F) +
    xlab(paste0(&quot;PC1: &quot;, round(percent_var[1] * 100), &quot;% variance&quot;)) + 
    ylab(paste0(&quot;PC2: &quot;, round(percent_var[2] * 100), &quot;% variance&quot;)) + 
  ggtitle(&quot;method: vegan::rda()&quot;)

p2</code></pre>
<p><img src="/post/2018-11-02-comparison-of-two-pca-methods_files/figure-html/pca2-1.png" width="768" /></p>
</div>
<div id="放在一起比较一下" class="section level1">
<h1>放在一起比较一下</h1>
<p>可以看出，虽然两种方法计算的数值有差异，但是坐标位置是一致。</p>
<pre class="r"><code>cowplot::plot_grid(p1,p2,labels = &quot;AUTO&quot;,ncol=1,align = &quot;hv&quot;)</code></pre>
<p><img src="/post/2018-11-02-comparison-of-two-pca-methods_files/figure-html/pca_plot-1.png" width="768" /></p>
</div>
